Modified detrended fluctuation analysis based on empirical mode decomposition 
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Detrended fluctuation analysis (DFA) is a simple but very efficient method for investigating 
the power-law long-term correlations of non-stationary time series, in which a detrending step is 
necessary to obtain the local fluctuations at different timescales. We propose to determine the local 
trends through empirical mode decomposition (EMD) and perform the detrending operation by 
removing the EMD-based local trends, which gives an EMD-based DFA method. Similarly, we also 
propose a modified multifractal DFA algorithm, called an EMD-based MFDFA. The performance of 
the EMD-based DFA and MFDFA methods is assessed with extensive numerical experiments based 
on fractional Brownian motion and multiplicative cascading process. We find that the EMD-based 
DFA method performs better than the classic DFA method in the determination of the Hurst index 
when the time series is strongly anticorrelated and the EMD-based MFDFA method outperforms 
the traditional MFDFA method when the moment order q of the detrended fluctuations is positive. 
We apply the EMD-based MFDFA to the one-minute data of Shanghai Stock Exchange Composite 
index, and the presence of multifractality is confirmed. 

PACS numbers: 05.45.Tp, 05.40.-a, 05.45.Df, 89.65.Gh 



I. INTRODUCTION 

The dynamics of an evolving complex system can be 
usually recorded as time series, whose temporal cor- 
relation structure embeds much information about the 
interactions among the microscopic constituents of the 
system. There are a wealth of approaches proposed 
to determine the correlation strength of a time series 
[l|, 0. In recent years, the detrended fluctuation anal- 
ysis (DFA) has become the most extensively adopted 
method. The idea of DFA was originally invented to in- 
vestigate the long-range dependence in coding and non- 
coding DNA nucleotide sequences 0, H|. An extension 
of the DFA method can be used to unveil the multifrac- 
tal nature hidden in time series, termed the multifractal 
DFA (MFDFA) @. The (MF)DFA method can also be 
generalized to investigate higher-dimensional fractal and 
multifractal measures @. 

There are also other applications of the DFA method. 
One example is the investigation of lagged correlations of 
nonstationary signals, in which the so-called lagged DFA 
can determine the largest correlation uncovering the ex- 
istence of underlying delays in the evolution of real time 
series 0. On the other hand, there are many situations 
where several variables are simultaneously recorded that 
exhibit long-range dependence or multifractal nature. A 



'Electronic address: xyqian@ecust.cdu.cn 
^Electronic address: wxzhou@ccust.cdu.cn 
t Electronic address: gfgu@ccust.cdu.cn 



detrended cross-correlation analysis (DXA) was proposed 
to investigate the long-range cross-correlations between 
two nonstationary time series, which is a generalization 
of the DFA method [l,©]. These methods can be easily 
generalized to study the power-law cross-correlations in 
higher-dimensional fractal or multifractal signals. 

A common step of all the aforementioned methods is 
to remove the local trends at different timescales. How- 
ever, for a real time series, it is usually not known if there 
is a trend and, if any, what the functional form of this 
trend is. A conventional strategy is to assume that the 
local trends are in the form of polynomial, which works 
quite well. Recently, a model-free timescale-adaptive de- 
trending approach has been proposed, which is based on 
the empirical mode decomposition (EMD) approach [l(| • 
For a given time series, EMD is designed to decompose it 
into a limited number of intrinsic mode functions (IMFs) 
and a residue component. The EMD-based trend can be 
considered as the combination of the monotonic residue 
and IMFs which are significantly distinguishable from the 
pure white noise, and the timescalc of an EMD-based 
trend is the averaged timescale of the decomposed IMFs. 
In this work, we employ the EMD as a detrending tool 
to modify the DFA algorithm, as well as the MFDFA 
algorithm. 

There are several relevant studies aiming at applying 
the EMD to study the correlations in time series. Nu- 
merical experiments based on fractional Gaussian noise 
with Hurst index Hq show that the EMD acts as a dyadic 
filter bank which is able to extract the Hurst index H in 
good agreement with the "true value" H when H ^ 0.5 
but shows a clear deviation from Hq when Hq < 0.5 
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A close scrutiny unveils that the extracted H is 
systematically greater than H [ll|. Another idea is to 
perform EMD on the original time series to remove the 
trend or seasonality and a DFA is followed, which was 
applied to study the correlation properties of long daily 
ozone records [12[. A slight different method is to con- 
duct the EMD dctrending process and the DFA directly 
on the original time series and a two-parameter scale of 
randomness for the DFA is proposed to replace the DFA 
scaling exponent, which was applied to identify charac- 
teristics and complexity of heartbeat interval caused by 
the effects of aging and of illness [l3|. Indeed, there is 
no power-law dependence of the detrended fluctuation 
function on the timescale. 

As we will show below, our EMD-bascd DFA and 
MFDFA algorithms arc different from the previous efforts 
in the combination of EMD and DFA. The performance 
of the EMD-based DFA and MFDFA methods is assessed 
with numerical experiments based on fractional Brown- 
ian motion and multiplicative cascading process. We find 
that the EMD-based DFA method performs better than 
the classic DFA method in the determination of the Hurst 
index when the time series is strongly anticorrelated, and 
the EMD-bascd MFDFA method outperforms the tradi- 
tional MFDFA method when the moment order q of the 
detrended fluctuations is positive. 

The paper is organized as follows. In Sec. [TH we re- 
view the algorithms of DFA, MFDFA and EMD and pro- 
pose the EMD-based DFA and MFDFA algorithms. In 
Sec. IIII1 we test the performance of the EMD-based DFA 
and MFDFA algorithms with extensive numerical exper- 
iments. The EMD-based MFDFA method is applied in 
Sec. IIVI to investigate the multifractal nature of the re- 
turn time series of the Shanghai Stock Exchange Com- 
posite index. Section [V] gives a brief summary. 

II. THE EMD-BASED DFA AND MFDFA 
ALGORITHMS 

A. The DFA algorithm 

The original DFA algorithm contains the following five 
steps |, || . 

Step 1. Consider a time series x(t), t = 1,2, ■■ • , N. 
First construct the cumulative sum 

t 

u(t) = 5>(<), t = l,2,---,N . (1) 

i=i 

Step 2. The new series u(t) is partitioned into iV s dis- 
joint segments of the same size s, where N s = [N/s]. 
Each segment can be denoted by u v such that u v (i) = 
u(l + i) for 1 ^ i ^ s, where I = (v — l)s. 

Step 3. In each segment u Vl we determine the local 
trend u v with the method of polynomial fitting. When 
a polynomial of order £ is adopted in this step, the DFA 
method is called DFA-^ (DFA-1 if £=l, DFA-2 if £=2, 



DFA-3 if £=3, and so on). We can then obtain the resid- 
ual sequence 

e v (i) = u v (i) - u v (i) , l^i^s. (2) 

Step 4- The detrended fluctuation function F(v, s) of 
the segment u v is defined as the root of the mean squares 
of the sample residuals e v (i) 

[F(v,s)} 2 = - ^2[e v (i)} 2 . (3) 
s z — ' 

i=l 

The overall detrended fluctuation is calculated by aver- 
aging over all the segments, that is, 

1 Ns 

[ F ( s )] 2 = ^rE[ F ( v ' s )] 2 ■ ( 4 ) 

S D — 1 

Step 5. Varying s, we can determine the power-law 
relation between the detrended fluctuation function F(s) 
and the timescale s, 

F(s) ~ s H , (5) 

where H is the DFA scaling exponent. In many cases in- 
cluding the fractional Brownian motions, the DFA scaling 
exponent H is identical to the Hurst index P, Q , which is 
related to the power spectrum exponent r\ by r\ = 2H — 1 
[Til , [l5| and thus to the autocorrelation exponent 7 by 
j = 2-2H. 

B. The MFDFA algorithm 

The MFDFA method has the same first three steps as 
the DFA method, and we need only to revise the last two 
steps @. 

Step 4- The gth order overall detrended fluctuation is 
calculated as follows, 

where q can take any real value except for q = 0. When 
q = 0, we have 

F ( S )=exp|^-^ln[F(v !S )]| , (7) 

according to L'Hospital's rule. 

Step 5. Varying the value of s, we can determine the 
power-law dependence of the detrended fluctuation func- 
tion F q (s) on the size scale s, which reads 

F q (s) ~ s h ^ , (8) 
where h(q) is the generalized Hurst index. 
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It is obvious that the DFA is a special case of the 
MFDFA when q = 2. In the standard multifractal for- 
malism based on partition function, the multifractal na- 
ture is characterized by a spectrum of scaling exponents 
r(g), which is a nonlinear function of q [l6|. For each q, 
we can obtain the corresponding traditional r(q) function 
through 

T(q) = qh(q) - D f , (9) 

where Df is the fractal dimension of the geometric sup- 
port of the multifractal measure. 



Note that, the trend r n (i) should be determined for each 
segment separately at each timescale. 

This gives the EMD-based DFA and MFDFA methods. 
One can see that these methods differ essentially from the 
ones proposed in the previous works [TH, [l3[ . 



III. VALIDATING THE METHODS THROUGH 
NUMERICAL EXPERIMENTS 

A. EMD-based DFA of fractional Brownian 
motions 



C. The EMD algorithm 

Empirical mode decomposition is an innovative data 
processing algorithm for nonlinear and non-stationary 
time series [lfj. It decomposes the time series x{t) into 
a finite number of intrinsic mode functions, which sat- 
isfy the following two conditions: (1) in the whole set 
of data, the numbers of local extrema and the numbers 
of zero crossings must be equal or differ by 1 at most; 
and (2) at any time point, the mean value of the "upper 
envelope" (defined by the local maxima) and the "lower 
envelope" (defined by the local minima) must be zero. 

The decomposing process is called a sifting process, 
which can be described with the following six steps [l(|: 
(1) Identify all extrema of x(t); (2) Interpolate the local 
maxima to form an upper envelope U(x); (3) Interpo- 
late the local minima to form a lower envelope L(x); (4) 
Calculate the mean envelope: 

H(t) = [U(x) + L(x)]/2; (10) 

(5) Extract the mean from the signal 

g(t) = x(t) - fi(t); (11) 

and (6) Check whether g(t) satisfies the IMF conditions. 
If YES, git) is an IMF, stop sifting; If NO, let x(t) = g(t) 
and keep sifting. Finally, we obtain 

n 

r n (t) = x(t)-^2 9i (t), (12) 
i=i 

where r„ is a residue representing the trend of the time 
series. 

D. The EMD-based DFA and MFDFA algorithms 

We can now embed the EMD algorithm into the DFA 
and MFDFA to modify the third step of the algorithms, 
while keeping all other steps unchanged. 

Step 3. For each segment u v , we obtain the EMD- 
based local trend u v = r n {i) with the sifting process. We 
can then obtain the residuals 

e v (i) = u v (i) — r n (i) , l^i^s. (13) 



We test the EMD-based DFA with synthetic fractional 
Brownian motions (FBMs). In this paper, we use the free 
MATLAB software FracLab 2.03 developed by INRIA to 
synthesize fractional Brownian motions with Hurst index 
Hq. In our test, we investigate fractional Brownian mo- 
tions with different Hurst indices H ranging from 0.1 to 
0.9 with an increment of 0.1. The size of each time series 
is 2 16 = 65536. For each Ho, we generate 100 FBM time 
series. Each time series is analyzed by the EMD-based 
DFA algorithm. 

0i . ■ ■ 1 
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FIG. 1: Log-log plots of the detrended fluctuation function 
F(s) with respect to the timescale s for two randomly selected 
FBM time series with H = 0.2 and H = 0.8 using the EMD- 
based DFA algorithm. The solid lines are the least-squares 
power-law fits to the data. 

In Fig. [TJ we show the log-log plot of the detrended 
fluctuation F(s) as a function of the timescale s for two 
randomly selected synthetic fractional Brownian motions 
with Hq = 0.2 (strongly anticorrelated) and Hq = 0.8 
(strongly correlated), respectively There is no doubt 
that the power-law scaling between F(s) and s is very ev- 
ident and sound, and the scaling range spans more than 
two orders of magnitude. The estimates of the Hurst in- 
dices are H = 0.197 and H = 0.796 for the anticorrelated 
and correlated time scries, which are very close to the cor- 
responding Hq values. The EMD-based DFA algorithm 
is able to well capture the self-similar (or self-affine) na- 
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ture of the fractional Brownian motions and results in 
precise estimation of the Hurst index. 

We confirm that there is also a power-law dependence 
of F(s) on s for other synthetic FBMs with different 
Hurst index Ho- For each Ho, we determine the Hurst 
index H for each FBM time series as done in Fig. [T] and 
obtain 100 H values. The mean of the 100 H values is 
calculated for each Ho- The resultant mean Hurst in- 
dices H are plotted against H in Fig. [2j We can see 
that the estimated Hurst indices H are very close to the 
preset values H . The deviation of the estimated Hurst 
index H from Ho becomes larger for larger values of Ho- 
When comparing with the EMD variance method which 
fails to give the estimates of Hurst index when Ho < 0.5 
[ll[ , we find that the EMD-based DFA method performs 
significantly better. 




FIG. 2: Assessing the performance of the EMD-based DFA 
method through extensive numerical experiments with frac- 
tional Brownian motions. A comparison of the estimated 
Hurst index H with the true value Ho is illustrated. The 
errorbar is determined by the standard deviation of the 100 
estimated H values for each Ho. 

Although Fig. [2] shows that the EMD-based DFA 
method is able to determine the Hurst index of a given 
fractional Brownian motion with high accuracy, it is nec- 
essary to compare the performance of the EMD-based 
DFA algorithm with the original DFA algorithm. We 
thus perform the DFA-1, DFA- 2 and DFA-3 algorithms 
on the same group of the FBM time series that were 
analyzed in Fig. [2] The relative ratio of the estimated 
Hurst index H with reference to the true value Hq is 
determined for each Hq for each of the four algorithms. 
Figure [3] digests the relative ratio H/Ho as a function 
of H for the four algorithms. For the EMD-based DFA 
algorithm, H/Ho increases smoothly and its maximum 
is less than 1.03. It means that this method overesti- 
mates H less than 3%. In contrast, the relative ratio 
H/Ho of the DFA-£ algorithms decreases with Ho and 
increases with the polynomial order £. For small Hq, the 
EMD-DFA method significantly outperforms the DFA-/? 



methods. For DFA-3, the relative deviation is as large as 
15%. For large Ho, the DFA-£ methods do a better job 
than the EMD-based DFA method. 
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FIG. 3: Assessing the performance of the EMD-based DFA 
method with reference to the original DFA methods using 
numerical experiments on fractional Brownian motions. The 
relative ratio H/Ho of the estimated Hurst index H over the 
true value Ho is plotted as a function of Ho for the EMD- 
based DFA and the original DFA algorithms. 



B. EMD-based MFDFA of multifractal signals 

We now turn to assess the performance of the EMD- 
based MFDFA algorithm with synthetic multifractal sig- 
nals. There are many methods proposed to generate 
multifractal signals, such as the multiplicative cascading 
method [13, E, Il9|. t he fractionally integrated singu- 
lar cascade method [2fJ, [2l|, [22| , the random W cascades 
method [22, 23], and so on. The multiplicative cascading 
process is widely used to model multifractal measures in 
many complex systems. As the simplest one, the p model 
was originally invented to simulate the energy-dissipation 
field in turbulent flows In this work, we adopt the 
textbook p model to generate multifractal signals. 

We start from a line and partition it into two segments 
of the same length and assign two given proportions of 
measure p\ = 0.3 and p-2 = 1 — Pi to them. Then each 
segment is divided into two smaller segments and the 
measure is redistributed in the same multiplicative way. 
This procedure is repeated for 16 times and at last we 
generate a multifractal signal of size 2 16 = 65536. If the 
multiplicative cascade process goes to infinity, the mass 
exponent function has an analytic expression as follows 



-{q)=-Hp\+pl)/\n2 . 



(14) 



The empirical mass exponent function of the constructed 
multifractal signal should be well approximated by 

E q . m 
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We perform the EMD-based MFDFA on the bino- 
mial measure and determine the empirical mass expo- 
nent function r(q), which is illustrated in Fig.[4j We also 
draw the theoretical line, Eq. (fTl| . in Fig. 2] for compar- 
ison. It is clear that the two curves overlap with each 
other. In addition, we also perform the MFDFA on the 
same multifractal signal and obtain the empirical r(q) 
function, which is plotted in Fig. 0] as well. We observe 
a marked discrepancy between the empirical r(q) curve 
extracted based on the classical MFDFA method and the 
theoretical curve when q 2. In addition, we find that 
the classical MFDFA method systematically underesti- 
mates the r(q) values when q 2, which is also observed 
in the case of two-dimensional multifractal surfaces 
This test shows that the EMD-based MFDFA method is 
able to extract the multifractal nature of signals more 
accurate than the classical MFDFA method at least in 
certain situations. 
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FIG. 4: Plots of r(q) extracted from the EMD-based MFDFA 
and the classical MFDFA as a function of q. The continuous 
line is the theoretical formula (11411, 



and the size of the data is 471202. For comparison, wc 
have also conducted the classical MFDFA with the poly- 
nomial order I = 1,2,3. The resulting r(q) functions for 
—5 ^ q ^ 5 are illustrated in Fig. [5] It is evident that 
the four algorithms give consistent results since the four 
r(q) curves overlap. 
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FIG. 5: (Color online) Multifractal analysis of the 1-min high- 
frequency returns of the SSEC index from 4 January 2000 to 
18 April 2008 using the EMD-based MFDFA and the classical 
MFDFA methods. 

There are two intriguing characteristic points in Fig. [5] 
When q — 2, r(q) = for all the four curves. It follows 
that the Hurst index of the SSEC returns is H = h(2) = 
[t(2) + l]/2 = 0.5, which is in agreement with the well- 
known fact that stock returns are uncorrelated. When 
g = 0, Fig. [5] gives r(q) = —1, which is in line with 
Eq. ([9]). In certain sense, these two points verify the 
correctness of the four algorithms. 



SUMMARY 



IV. APPLICATION TO THE SSEC INDEX 

In this section, we apply the EMD-based MFDFA 
method to study the multifractal nature of the high- 
frequency return time series of the Shanghai Stock Ex- 
change Composite (SSEC) index. The multifractal prop- 
erties of financial returns have been investigated exten- 
sively [24l I25] ]. Concerning the Chinese stock market, 
there are also numerous multifractal analyses based on 
the multiplier method the MFDFA method [p,[2§, 
and the partition function approach [27], [H, l3(il. 1 3 ll 1321 
[33} ] . The presence of multifractality in the Chinese stock 
market is well documented. 

We have performed the EMD-based MFDFA on the 1- 
min high-frequency data of the SSEC index from 4 Jan- 
uary 2000 to 18 April 2008. There are 240 minutes in the 
double continuous auction on each trading days [34], |35| , 



In summary, we have proposed a modified detrended 
fluctuation analysis for fractal and multifractal signals 
based on the empirical mode composition. The polyno- 
mial local trend in the classic DFA algorithm is replaced 
by an EMD-based local trend. The modified (MF)DFA 
is called the EMD-based (MF)DFA. 

The performance of the EMD-based DFA and MFDFA 
methods is assessed with extensive numerical experi- 
ments based on fractional Brownian motion and multi- 
plicative cascading process. For the EMD-based DFA 
method, we investigated different Hurst index 0.1 ^ 
Hq ^ 0.9 and generated 100 fractional Brownian mo- 
tions for each Ho- The accuracy of the estimated Hurst 
indices H obtained from the EMD-based DFA with ref- 
erence to Hq is compared with that from the classical 
DFA. We found that the EMD-based DFA performs bet- 
ter than the classic DFA method in the determination 
of the Hurst index when the time series is strongly anti- 
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correlated, especially when Hq < 0.3, while the classical 
DFA outperforms when H > 0.4. In all cases, the EMD- 
based DFA is able to determine the Hurst index H with 
a deviation less than 3% from the true value Hq, that is 
(H - H Q )/H < 3%. 

For the EMD-based MFDFA method, we constructed 
a multifractal signal based on the p model. The empirical 
mass exponent functions r(q) of the EMD-based MFDFA 
and classical MFDFA were compared with the analyti- 
cal expression. We found that the EMD-based MFDFA 
outperforms the traditional MFDFA method when the 
moment order q of the detrended fluctuations is positive. 
The usefulness of the EMD-based MFDFA in the multi- 
fractal analysis is thus validated. 

As an empirical example, we have applied the EMD- 
based MFDFA to the 1-min return data of Shang- 
hai Stock Exchange Composite index. The EMD- 
based MFDFA gives very similar results as the classical 
MFDFA methods with different detrending polynomials. 
The presence of multifractality is confirmed. 



We conclude that the EMD-based DFA and MFDFA 
methods have comparable performance as the classical 
DFA and MFDFA methods in the analysis of fractal and 
multifractal time series. In certain cases, the EMD-based 
methods can give better results. The only shortcom- 
ing of the EMD-based DFA and MFDFA algorithms is 
that the detrending process based on EMD is more time- 
consuming. 
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